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A mathematical model of interacting species filling ecological niches left by the 
extinction of others is introduced. Species organize themselves into genera of all 
sizes. The size of a genus on average grows linearly with its age, confirming a general 
relation between Age and Area proposed by Willis. The ecology exhibits punctuated 
equilibrium. Analytic and numerical results show that the probability distribution 
of genera sizes, genera lifetimes, and extinction event sizes are the same power law 
P(x) ~ 1/x 2 , consistent with paleontological data. 
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Many years ago, Willis noted that genera could be composed of many species or of only 
one; he noted that similar regularities in the statistical properties of genera occured whether 
one is studying flowering plants or e.g. beetles. Attempting to formalize these observed 
regularites, he postulated a relation between the age and size of a genus, "Age and Area," 
which states that older genera on average include more species than younger ones . With 
Yule he noted a power law relation for the number of genera with s species, P gen (s) ~ s~ T 
with r approximately 2 Recently, Burlando |§ observed scaling behavior across the 

taxonomic hierarchy also giving r ~ 2 (see in addition Ref. 0]). Similarly, the distribution of 
life times, t , of fossil genera || can be described by a power law Pu f e (t) ~ 1 / t Tt with — 2 || . 
When viewed at sufficiently large time scales, the pace of extinction itself is episodic with long 
periods of stasis interrupted by sudden bursts of mass extinction 0. Punctuated equilibrium 
with scale-free extinctions has been attributed || to the self-organized critical |10[ dynamics 



of strongly interacting species, without the need for catastrophic exogenous causes such as 
meteorites. The punctuated equilibrium process governing large scale evolution may be 
sufficiently robust or universal to be captured by an abstract, mathematical model. Bak 
and Sneppen |TT| introduced a self-organized critical model of coevolving "species" where the 



least fit undergoes pseudo-extinction and affects the fitness of other species in the ecology, 
leading to extinction events of all sizes. However, their dimension-independent lifetime 
distribution P H f e (t) ~ 1/t is in disagreement with Raup's [|,|J paleontological data. 
More importantly, their model lacks any emergent taxonomic structure, which is an essential 
part of large scale organization in evolution. 

Here we show that these three, seemingly unrelated distributions for extinction event 
sizes, genera sizes, and genera lifetimes, which characterize the large scale behavior of evo- 
lution, can be unified in terms of a simple mathematical process. We introduce an abstract 
model for large scale evolution and study it both analytically and numerically. In our model, 
surviving species can diversify into ecological niches left by previous extinction. This is im- 
plemented in terms of a Polya urn type of process. All species are subject to a general drift 
over time to lower viability which eventually leads to their extinction. This is consistent 



with the view that most changes in the ecology have a deleterious effect on currently existing 
species. In addition, due to interactions between species they experience changes in their 
viability, either favorable or unfavorable, arising from previous extinction of other species in 
the ecology. Species in our model organize themselves into genera of all sizes. The emergent 
genera obey a general Age and Area |IJ relation which we find is linear (s ~ t) in contrast 
to Yule's conjecture Ins ~ t. This large scale organization of species provides a simple 
mechanism which shapes all three probability distributions into a power law P(x) ~ 1/x 2 , 
where x is the extinction event size, genera lifetime, or genera size. In all three cases the 
1/x 2 behavior is consistent with previously reported paleontological data f|Pf||3| , indicat- 
ing that our model may plausibly describe a universality class sufficiently broad to include 
real evolution. It can be tested further by directly comparing the age and size of extinct 
genera with our result that on average s ~ t. 

Our model was inspired, in part, by considering a more complicated "connection" model 
introduced by Sole and Manrubia |I2 |. The advantage of our model is that it is extremely 



simple and robust. Its simplicity makes it easier to study large systems numerically; it 
is also analytically tractable. That such a simple model exists which describes a process 
giving large scale organization in evolution together with the above mentioned distributions 
is significant we think because it illustrates a potentially universal mechanism that would 
apply even beyond the context of biological evolution discussed here. 

We begin by briefly describing the connection model. An N x N interaction matrix W 
defines the interaction, either favorable or unfavorable, between N objects that represent 
species. For a species i, the output elements Wij define its affect on the other species j, while 
its viability is the sum of its input elements v{i) = J2j Wji. If Wji > Wki then species j has a 
more beneficial effect (and species k has a more deleterious effect) on the ability of species i to 
survive. If the viability v (i) < 0, then species i goes extinct, and the connection elements of 
the rows and columns for that extinct species are replaced with a copy of the corresponding 
elements of another surviving species. This copying in turn changes the viability of other 
species, and leads to a chain reaction of extinction events. The system is driven by slow 



random changes in the matrix W which tend to lower the viability of the surviving species, 
and slowly differentiate copies from each other leading to speciation. They observed that 
the connection model exhibits extinction events of all sizes where the species that go extinct 
together tend to be recent copies. 

In our one-dimensional model we assign to each of N particles that represent species an 
integer viability v(i). The dynamics consists of three steps as illustrated in Fig. 1: (1) species 
drift stochastically to lower viability; (2) species with viability below a threshold v c become 
extinct. The extinct species are each replaced with a "daughter" speciation of a surviving 
species. This is the Polya urn mechanism in our model; (3) Due to interactions between 
species the surviving species receive a change in their viability resulting from the extinction 
event. Specifically, at each time step the following operations are performed in parallel for 
all species (i): (1) with probability 1/2, v(i) = v(i) — 1; otherwise v(i) is unchanged; (2) 
for each % such that v{%) < v c a surviving species (J) with v(j) > v c is selected at random 
and v(i) = v(j). This step represents a speciation event where one species branches into 
two. (3) all N — s species that survived extinction receive a coherent influence q(s), so that 
v (j) = v {j) + <?( s )- After an extinction event of size s, q(s) is chosen from the uniform 
distribution — s < q(s) < s. Thus, only large extinctions can cause large subsequent changes 
in the ecology. The form of q(s) is elaborated on later. It is important to note that, unlike 
the connection model, our model's behavior is robust with respect to varying the parameter 
v c , since the entire system is translationally invariant in viability. 

Our model can be viewed as an example of transport in one dimension, where particles 
are conserved. In the steady state the smooth drift of species toward lower viability will 
be balanced by the intermittent replacement of extinct species with speciations of surviving 
ones, which by definition have higher viability. The average viability in the system v = 
jfYU v {^) exhibits stick-slip behavior as shown in Fig. 2, similar to the behavior observed 
in the connection model ||12|| . 

Due to replacement of extinct species with speciations of surviving ones, the species tend 
to form groups with similar viability, which drift and diffuse together toward the extinction 
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threshold. The state of the system may be characterized by n(v,t), the number of species 
of viability v at time t, where = N. A snapshot of the system is shown in Fig. 

3. At a microscopic scale, n(v,t) is peaked with well defined bumps that give rise to the 
temporally intermittent sequence of extinctions as shown in the insert of Fig. 2. 

We can identify all species within each bump as members of the same genus for the 
following reasons: Each viability bump is separated from the others by an empty interval 
where n(v,t) = 0. Since these empty intervals cannot be filled by the replacement of 
extinct species with speciations of surviving ones, the dynamics tends to maintain the sharp 
separation between different bumps. Therefore, they are long-lived metastable entities. By 
making a histogram of the genera sizes, or area under each bump, observed in snapshots at 
spaced time intervals we find a power law for the number of species within each genus as 
shown in the insert of Fig. 3 with an exponent r ~ 2. Also, the total number of genera in 
the system displays an intermittent pattern of diversification (increase) and contraction in 
time, qualitatively similar to real data |p~3f . 

The age of a newly created species following extinction is set to zero, and incremented 
by one unit at each step in the simulation. The age of a genus is the age of the oldest species 
in the corresponding bump. When a bump passes through the extinction threshold v c , we 
measure its age t and size s (or area). The distribution of sizes of extinct genera is the same, 
within numerical accuracy, as the snapshot distribution described above. We numerically 
determined the relation between the age and size of extinct genera in a system of size 
N = 1000 including 10 7 time steps, and found a linear relation t = ms, with m ~ 0.6. This 
numerical result indicates that emergent genera on average grow at constant rate irrespective 
of their size. The numerical result r = 2, then implies r t = 2, in agreement with real 
data. Data collapse of the distribution of extinction event sizes for different system sizes 
also indicate a power law with exponent r ext = 2, as shown in Fig. 5. Note that at the 
beginning of the numerical simulation, with random initial conditions, there are only small 
extinctions and small genera. Thus the power law distributions observed in the steady state 
are "emergent" ; they are consequences of the self-organized critical dynamics of our model. 
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The first step in our model, drift to lower viability, takes into account slow random 
mutations of the matrix elements Wij in the connection model that tend to lower the via- 
bility of all species. This slow external driving, similar to that used in earthquake models 
||14j| , represents the cumulative effect on species of small changes in the environment, which 
we propose tend to make species less able to maintain their population over time. The 
second step represents true extinctions of species. The third step represents the effect of 
these extinctions on the surviving species. For more details see Ref. [fL5 l. Our preliminary 



results indicate that the connection model also exhibits emergent genera with a broad size 
distribution, intermittent diversification, as well as an Age and Area relation. 

We now discuss the analytic results for the transport model. In the stationary state, 
Pext{s) is the probability distribution to have an extinction event (avalanche) of size s and 
G(q) is the probability distribution to have an influence of size q. These distributions are 
self- consistently related via 

<?(?)= E §^4 (1) 

all s>q 

Pe^) = E G (^(E^) ) S ) , (2) 
all q v=0 

where n(v ) is the time averaged viability profile in the steady state. The first equation is 
exact. The second assumes that the avalanche distribution comes from influences, q, on the 
time average viability profile, rather than the actual time dependent profile. The extinctions 
and influences are treated in terms of their full probability distributions, while the viability 
profile of species is treated only in terms of its average. This can be justified a posteriori in 
terms of a separation of scales argument similar to singular diffusion |16| . 

Next, we assume that the cumulant of n(v ) is not singular around v = 0, so that it has 
a Taylor series expansion. In the interval 1 <C q <C N, where N — > oo, J Q n(v)dv = Aq + .... 
Combining Eqs. (1,2) with the Taylor series expansion gives 



dG{q) 



dq 



G(jj) for l«g«JV . (3) 



2Aq M 



It is easy to show that Eq. (3) has a scaling solution G(q) ~ q Text where ln(2r ext 



{Text — 1) In A. Also the avalanche and influence distributions are asymptotically the same; 
PeM) ~ G(q). 

The steady state equation for the time average profile is 



-i N-l s p l \ i N-l p l \ 

n{v) = EE + ?) + n{v + g + 1)) + n(«) + fi( w + 1)) E 



For large AT, we try the solution n(v) = n a e CV ^ N and find to leading order in N 



(4) 



1 N-l s p i \ N-l p l \ 

1 = ~ E E ^e"^(2 - c/iV) + E ■ (5) 

Expanding for q << N, the g = part of the first term on the right hand side gives 1 — jjy. 
The leading part cancels the number one on the left hand side of Eq. [|, and the negative 
remainder which comes from the drift must be cancelled by the remaining terms in the 
equation. Completing the expansion in q, only the even terms survive the symmetric sum 
over q. These terms are all positive as is the last term in Eq. [5| From Eq. (3) all of these 
positive terms scale ~ jv 1-7 "^' . Only when the positive terms cancel the only 

negative term (^)- In this case, a consistent solution exists for the exponential profile. Our 
numerical simulation results show that the average profile is indeed exponential with c ~ 6 
with A ~ 4, both confirming r ex t = 2. 

Note that our theory thus far has removed genera bumps by only treating the time 
average profile. The weak In iV divergence of both the average size of influences and average 
size of extinction events justifies our separation of scales assumption for large N. Since the 
shifts are small relative to N, only a finite number of genera on average pass the extinction 
threshold following an extinction. Then r = r ext = 2, in agreement with the numerical 
simulation result. Finally, previous interpretations of available data from the fossil record 
for extinction size distributions, genera (and higher order taxa) abundance distributions, 
and lifetime distribution of genera are consistent with our unified result in terms of a simple 
model that they are each decaying power laws with exponent 2. 
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FIGURES 

FIG. 1. Dynamics of the model. The horizontal axis is the viability and the blocks repre- 
sent species. The dotted line is the extinction threshold. (0) Initial configuration. (1) Leftward 
stochastic drift. (2) Extinction and replacement. (3) Coherent influence to the survivors. Here the 
extinction had size s = 5 and the influence of the extinction had value q = —2. The species that 
move at each step are shaded. 

FIG. 2. The average viability as a function of time in a system of size N = 1000 exhibiting 
stick-slip dynamics. The steep jumps, or slip events, are followed by slow relaxation to the threshold 
for extinction. The insert shows the temporal sequence of extinction event sizes over the same 
interval. 

FIG. 3. A snapshot of the viability profile n(v, t) in a system of size N = 400. The pattern 
is intermittent with both small and large bumps. The insert shows the distribution of genera sizes 
averaged over a total time interval of 10 7 steps with a snapshot taken every 100 time steps for a 
system of size N = 1000. The curve can be described as power law with a cutoff at the system 
size. 

FIG. 4. Data collapse result. P(s, N) is the probability to have an extinction event of size 
s in a system of size N. The plateaus for different system sizes show that P(s,N) = F(s/N)/s 2 
where F(x) is simple scaling function which is constant for x <C 1 and approaches zero as x — ► 1. 
This agrees with our analytic result that r ex t = 2. 
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